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^ ■ Abstract 

■ , We explore the utility of the recently proposed alpha equations in provid- 

ed , 
c~| , ing a subgrid model for fluid turbulence. Our principal results are comparisons 
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of direct numerical simulations of fluid turbulence using several values of the 
parameter alpha, including the limiting case where the Navier-Stokes equa- 
tions are recovered. Our studies show that the large scale features, including 
statistics and structures, are preserved by the alpha models, even at coarser 
resolutions where the fine scales are not fully resolved. We also describe the 
differences that appear in simulations. We provide a summary of the principal 
features of the alpha equations, and offer some explanation of the effective- 
ness of these equations used as a subgrid model for three-dimensional fluid 
turbulence. 



I. INTRODUCTION 

Holm et al. JT|, introduced the "alpha models" for the mean motion of ideal incom- 
pressible fluids as the n-dimensional generalization of the one- dimensional Camassa-Holm 
(CH) equation. The ID CH equation describes shallow water waves with nonlinear dispersion 
and admits soliton solutions called "peakons" ||. Its n-dimensional generalization describes 
the slow time dynamics of fluids in which nonlinear dispersion accounts for the effects of the 
small scale rapid variability upon the mean motion. The fluid transport velocity is found 
by inversion of a Helmholtz operator acting on fluid circulation (or momentum) velocity. 
This operator contains a length scale that corresponds to the magnitude of the fluctuation 
covariance; the application of this operator smoothes the transport velocity relative to the 
circulation velocity. This length scale is denoted by a in references |jl|, ||, ||, hence the 
name alpha models for these mean fluid motion theories. The alpha models for self consistent 
mean fluid dynamics are derived by applying temporal averaging procedures to Hamilton's 
principle for an ideal incompressible fluid flow. The resulting mean fluid motion equations 
are obtained by using the Euler-Poincare variational framework Q, [0. (Euler-Poincare 
equations are the Lagrangian version of Lie-Poisson Hamiltonian systems.) Therefore, these 
equations possess conservation laws for energy and momentum, as well as a Kelvin-Noether 
circulation theorem that establishes how the time average (or perhaps statistical) properties 
of the fluctuations affect the circulation of the mean flow. These ideal fluid equations also 
describe geodesic motion on the volume-preserving diffeomorphism group for a metric con- 
taining the H 1 norm of the mean fluid velocity. Their geometrical properties are discussed 
in Q. Their relation to Eulerian and Lagrangian mean fluid theories is discussed in ||. In 
recognition of their origins, these mean fluid motion equations may be known equally well 
by either the name alpha models, or CH equations. 

Chen et al. |J [|J introduced phenomenological viscosity into the CH equation and 
proposed the resulting viscous Camassa-Holm equation (VCHE) as a closure approxima- 
tion for the Reynolds averaged equations of the incompressible Navier-Stokes fluid. They 



tested this approximation on turbulent channel and pipe flows with steady mean, by find- 
ing analytical solutions of the VCHE for the mean velocity and the Reynolds shear stress 
and comparing them with experiments ||. They found that the steady VCHE profiles are 
consistent with data obtained from mean flow turbulence measurements in most of the flow 
region for channels and pipes at moderate to high Reynolds numbers. Thus, Chen et al. 
demonstrated a connection between turbulence and the VCHE for steady, or mean solu- 
tions. In fact, the time- dependent VCHE in a periodic box has unique classical solutions 
and a global attractor whose fractal dimension is finite and scales according to Kolmogorov's 
estimate, N ~ (L/£d) 3 , where id = (z/ 3 /e) 1//4 is the Kolmogorov dissipation length [10|. We 



note that the time-dependent VCHE model is not equivalent to the Navier-Stokes equations 
with hyperviscosity, see @ -- ||, [K|- The VCHE model is also known as the Navier-Stokes 
alpha model, or the viscous alpha model. 

The purpose of this paper is to compare the statistics and structures of the velocity 
and vorticity fields at moderate Reynolds numbers in a direct numerical simulation (DNS) 
of the viscous alpha model with the corresponding results for the Navier-Stokes equations. 
Through this comparison, we hope to determine whether the Navier-Stokes alpha model 
can be used as a subgrid model for fluid turbulence. The paper is arranged as follows: the 
mathematical background of the viscous alpha model and its connections with the Navier- 
Stokes equations are given in the next section. The remainder of the paper compares the 
classical Navier-Stokes DNS turbulence results with those for the viscous alpha models. In 
Section III we present the numerical methodology and energy spectra. Stretching dynamics, 
including the vorticity structure and the alignment phenomena are presented in Section 
IV. In Section V we present various two-point statistics, including the probability density 
function (PDF) of the velocity increment and the velocity structure functions. Concluding 
remarks are given in Section VI. 



II. EULERIAN MEAN THEORY AND ALPHA MODEL EQUATIONS 

Holm, Marsden, and Ratiu [|l| , || used variational asymptotics to obtain evolution equa- 
tions for the Eulerian mean hydrodynamic motion of ideal incompressible fluids, employing 
an approximation of Hamilton's principle for Euler's equation in a Euclidean space setting. 
The method assumes that the Euler flow may be decomposed into its mean and fluctuating 
components at a fixed position in space. In their approach, a first-order Taylor expansion 
in the fluctuation amplitude is used to approximate the velocity field with the result that 
the L 2 metric in the Hamilton's principle giving rise to the Euler equations is replaced by 
an H l metric that produces the evolution equation of the Eulerian mean flow. We shall call 
this evolution equation the Euler alpha model, or the n-dimensional CH equation. 

Holm et al. [|j] give a geometrically intrinsic (i.e., coordinate free) derivation of these 
averaged equations by the procedure of variational asymptotics, namely, by deriving an 
averaged Lagrangian and using this Lagrangian to generate the equations via Hamilton's 
principle. This intrinsic setting is useful because many interesting flows, e.g., flows on spheres 
such as those in geophysics, do occur on manifolds. The standard decomposition into mean 
and fluctuating components is an additive decomposition, only valid in the presence of a 
vector space structure. We follow Holm || in presenting the derivation in Euclidean space. 

We develop the Euler-Poincare theory of advected fluctuations from the viewpoint of 
Eulerian averaging. Our point of departure is a Lagrangian comprised of the fluid kinetic 
energy in the Eulerian description, in which volume preservation is imposed by a Lagrange 
multiplier P (the pressure), 

LH = y t / 3 x|y|U(x,t;^)| 2 + P(x,t;^)(l- J D(x,t;^))| , (1) 

where D is the Eulerian volume element. We assume that there are two time scales for the 
motion: the fast time u and the slow time t. The traditional Reynolds decomposition of 
fluid velocity into its fast and slow components is expressed at a given position x in terms 



of the Eulerian mean fluid velocity u as 

U(x,*;w) = u(x,t)+u'(x,t;a;). (2) 

Following Holm [[J, we assume the Eulerian velocity fluctuation u'(x, t;u>) is related to an 
Eulerian fluid parcel displacement fluctuation — denoted as £(x, t;u>) --by 

dC 

^ + u-VC = C-Vu + u'(x,i; W ). (3) 

For purely Eulerian velocity fluctuations as in Eq. (0), this relation separates into two 

relations: the "Taylor-like" hypothesis of ||, 

dC 

^ + u-VC = 0; (4) 

and the relation ||, [|] 

= C- Vu + u'(x,t;w). (5) 

Hence, the Reynolds velocity decomposition (|j) separates the Lagrangian ([[]) into its mean 
and fluctuating pieces as 

L(u)= /'rf 3 x|^|u(x,t) + u , (x,t;cu)| 2 + P(x,t)(l- J D(x,t)U . (6) 

No modification is needed in the pressure constraint in this Lagrangian, because the Eulerian 
mean preserves the condition that the velocity be divergence free; hence, V-u = 0. It 
remains only to take the Eulerian mean of this Lagrangian, in which we assume (£ ) = 0. 
The Eulerian mean averaging process at fixed position x is denoted ( • ) with, e.g., 

i r T 

u(x,t) = (U(x,t;w)) = lim - / U(x, t ; w) duj . (7) 

T^oo T J 

By Eq. (^), the Eulerian mean kinetic energy due to the velocity fluctuation satisfies 

(|u'| 2 > = (C fc C'>u, fc • u, . (8) 

Thus, we find the following Eulerian mean Lagrangian, 

(L) = Jd 3 x | ^ [| u (x,t)| 2 + (C fe C')u,fe • u,i] + ^(x,t)(l - D(x,t)) } . (9) 



The advection relation (f|) implies the same advective velocity for each component of the 
symmetric Eulerian mean covariance tensor (( k ( l ). Thus, we have 

(| + u -v)(cV) = o. (10) 

Together, this relation and the continuity equation for the volume element D, 

3D 

^ r + V.Du = 0, (11) 

complete the auxiliary equations needed for deriving the equation of motion for the Eulerian 
mean velocity u from the averaged Lagrangian (L) in (0) by using the Euler-Poincare theory. 
The results of JT] allow one to compute the Euler-Poincare equation for the Lagrangian 
(L) in (|6[) depending on the Eulerian mean velocity u, and the advected quantities D and 
(C fc C) as 

( d ,9\1 SIL) 1 SIL) ,- 

_il, 1 S(L) d kl 
dx i 5D D 5(( k ( l ) dx l ^ ^ h 

We compute the following variational derivatives of the averaged approximate Lagrangian 

(L) in Eq. © 



^ = -5(* D <««)"^ 






( u , fe • u,j) . (13) 



${( k ( l ) 2 
The Euler-Poincare Eq. (|12|) for this averaged Lagrangian takes the form, 

-^ + u -V v + BjVtt* + VP W = - - (u fc • u,,)V (C fc C') , (14) 

where v = u-A B u with A D = ^-(d k D (( k ( l )di) and V-u = 0. (15) 

Its definition as a variational derivative indicates that v is a specific momentum in a certain 
sense dual to the velocity u. For more discussion of physical interpretations of u and 
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v, see ||. The Euler-Poincare equations fll4|) - ( fLq) define the Eulerian mean motion 
(EMM) model. Incompressibility of the Eulerian mean velocity u follows from the continuity 



equation (II ) and the constraint 5{L)/5P = 0. A natural set of boundary conditions is 

v-n = , u = , and n • (££) =0, on a fixed boundary. (16) 

Then, provided the Helmholtz operator 1 — Ao for D — 1 may be inverted, the Eulerian 
mean pressure P may be obtained by solving an elliptic equation. 

A. Reducing the EMM equation to the n-dimensional CH equation 

When the Eulerian mean covariance is isotropic and homogeneous, so that (C fc C0 = ot 2 8 kl 
(for a constant length scale a, whose magnitude is set by the initial conditions for the 
Eulerian mean covariance) then the EMM equation (|H|) reduces to the n-dimensional 
Camassa-Holm equation introduced in [|TJ, 0, namely, 

<9v 

- + u-Vv + «,W+VP w =0, V-u = = V-v, (17) 

at 

where v = u-ct 2 Au, P tot = P - -|u | 2 - — |Vu | 2 . (18) 

This n-dimensional CH equation set is an invariant subsystem of the Euler-Poincare system 
([14]), with definition (|l^) and advection law fllOD , because the homogeneous isotropic initial 
condition (C fc O = & 2 5 kl is invariant under the dynamics of equation (l0|). Hence, any of the 
formulae above remain valid if we set (( k ( l ) = a 2 5 kl , with constant a. 



Equations of the type ( |TTD but with additional dissipative terms were considered previ- 



ously in the theory of second grade fluids [1 1] and were treated recently in the mathematical 



literature fl2f , ||13|| . Second grade fluid models are derived from continuum mechanical prin- 



ciples of objectivity and material frame indifference, after which thermodynamic principles 
such as the Clausius-Duhem relation and stability of stationary equilibrium states are im- 
posed to restrict the allowed values of the parameters in these models. In contrast, the CH 
equation ( |T7| ) is derived here by applying asymptotic expansions, Eulerian means, and an as- 
sumption of isotropy of fluctuations in Hamilton's principle for an ideal incompressible fluid. 



This derivation provides the interpretation of the length scale a as the typical amplitude of 
the rapid fluctuations whose Eulerian mean is taken in Hamilton's principle. 

The n-dimensional CH equation ( |TTD implies the conservation of energy | J d 3 xu • v 
and helicity | J d 3 x v • curlv. Its steady vortical flows include the analogs of the Beltrami 
flows curlv = Au. In the periodic case, we define Vk as the k-th Fourier mode of the specific 
momentum v = (1 — a 2 A)u; so that v^ = (1 + a 2 |k| 2 )uk. Then Eq. (|17|) becomes flU, |2| 



n± — v k - i V Vp x (n x v n ) 

\ at ^-^ 1 + crpr / 

V p+n=k '^' / 



0, (19) 

— ± -r u!"iur / 

p+r 

where n^ is the Leray projection onto Fourier modes transverse to k (this ensures in- 
compressibility) . Hence, the nonlinear coupling among the modes is suppressed by the 
denominator when l + a 2 |p| 2 S>|n|. 

An essential feature of the n-dimensional CH equation (|TT| ) is that its specific momentum 
v is transported by a velocity u that is smoothed, or filtered, by application of the inverse 
elliptic Helmholtz operator (1 — a 2 A). The effect on length scales smaller than a is that 
steep gradients of the specific momentum v tend not to steepen much further, and that 
thin vortex tubes tend not to get much thinner as they are transported. Furthermore, as 
our present numerical simulations shall verify, the effect on length scales larger than a is 
negligible. Hence, the n-dimensional CH equation preserves the assumptions under which it 
is derived. 

B. Physical interpretation of v as the Lagrangian mean velocity 
The Stokes mean drift velocity is defined by fT4|| , 



(U) 5 =(C-V U '>. (20) 

Hence, Eq. (|5|) implies 

(U) 5 = - (c • vc • V)u = - Au + o(|C| 2 ) , (21) 



where 

A=(d,(C fc CW=A D | D=1 , (22) 

and we argue that V • C — °(IC| 2 )- Thus, we find that v satisfies, to order o(|£| 2 ), 

v = u-Au = u+(U) s = (U) L . (23) 

Therefore to this order, v in the EMM theory is the Lagrangian mean velocity. 

C. Kelvin circulation theorem for EMM and CH equations 



Being Euler-Poincare, the Eulerian mean motion (EMM) equation (|T4]) and its invariant 
reduced form the CH equation ([17]) for (C fc C') = o?b kl has a corresponding Kelvin-Noether 
circulation theorem, 

*l v.dx=-i/"/ V(u fc .Uj)xV<C fc C'MS, (24) 

at Jj(u) Z J JS(u) 

for any closed curve 7(11) that moves with the Eulerian mean fluid velocity u and surface 
S(u) with boundary 7(11). Thus in this Kelvin-Noether circulation theorem, the presence 
of spatial gradients in the Eulerian mean fluctuation covariance (C fc C') creates circulation of 
the Lagrangian mean velocity v = u — Au. 

D. Vortex stretching equation for the Eulerian mean model 

In three dimensions, the EMM equation ( |14D may be expressed in its equivalent "curl" 
form, as 

— v - u x (V x v) + V(P tot + u • v) = - - (u fc • u,,) V(C fc C'> , V • u = . (25) 

The curl of this equation in turn yields an equation for transport and creation for the 
Lagrangian mean vorticity, q = curlv, 

^ + u- Vq = q- Vu- -V(u fc -u/) x V(C fc C'), where q = curlv, (26) 



and we have used incompressibility of u. Thus u is the transport velocity for the gener- 
alized vorticity q, and the expected vortex stretching term q -V u is accompanied by an 
additional vortex creation term proportional to the Eulerian mean covariance gradient. Of 
course, this additional term is also responsible for the creation of circulation of v in the 
Kelvin- Noether circulation theorem ( ffij|) and vanishes when the Eulerian mean covariance is 
homogeneous in space, thereby recovering the corresponding result for the three dimensional 
CH equation Q, @. 

E. Energetics of the Eulerian mean model 

Noether's theorem guarantees conservation of energy for the Euler-Poincare equations 
(fUj), since the Eulerian mean Lagrangian (L) in Eq. (^ has no explicit dependence on time. 
This constant energy is given by 



E, 



i J d 3 x(\u | 2 + (CV)u, k • u, ; ) =^fd 3 xu-v. (27) 



Thus, the total kinetic energy is the integrated product of the Eulerian mean and Lagrangian 
mean velocities. In this kinetic energy, the Eulerian mean covariance of the fluctuations 
couples to the gradients of the Eulerian mean velocity. So there is a cost in kinetic energy 
for the system either to increase these gradients, or to increase the Eulerian mean covariance. 

F. Momentum conservation — stress tensor formulation 

Noether's theorem also guarantees conservation of momentum for the Euler-Poincare 
equation ([14]), since the Eulerian mean Lagrangian (L) in Eq. (j6|) has no explicit spatial 
dependence. In momentum conservation form, Eq. ( |14D becomes 

^---^-{y^+P8i-u k .nAC k C j }). (28) 

The boundary conditions are given in Eq. (fig 1 ). 
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G. A second moment turbulence closure model for EMM 



When dissipation and forcing are added to the EMM motion equation fl"H|) by using 



the phenomeno logical viscosity z/Av and forcing F, one finds a second moment Eulerian 
mean turbulence model given by 

(- + u -V )v + Vj Vv? + VP tot + - (u fe • u,,) V(CV> (29) 

= v Av + F , where V • u = , 

with viscous boundary conditions v = 0, u = at a fixed boundary. Note that the Eulerian 
mean fluctuation covariance (C fe C J ) appears in the dissipation operator A. In the absence 
of the forcing F, this viscous EMM turbulence model dissipates the energy E in Eq. Q27D 
according to 



— - = -v I d 6 x 
dt 



I d 3 x tr(Vu T • (CC) -Vu) + Au • Au . (30) 



This negative definite energy dissipation law is a consequence of adding viscosity with A, 
instead of using the ordinary Laplacian operator. In the isotropic homogeneous case of this 
model, where (( k ( l ) = a 2 5 M (for a constant length scale a) we find the viscous Camassa- 
Holm equation (VCHE), or the Navier-Stokes alpha model, 

(— + u •v]v + w i Va i + VP tot = ra 2 Av + F, where V-u = 0, (31) 

where A is the usual Laplacian operator for this case and v and P tot are defined in fll8D . 



H. Constitutive interpretation of the VCHE 

Chen et al. 0- M gave a continuum mechanical interpretation to the VCHE closure 



model, by rewriting the VCHE ([H]) in the equivalent constitutive form, 

-^ = divT , T = -pi + 2i/(l - a 2 A)D + 2a 2 D , (32) 

with V-u = 0,D = (l/2)(Vu+Vu T ), ft = (l/2)(Vu- Vu T ), and co-rotational (Jaumann) 
derivative given by D = dD/dt + DO — f2D, with d/dt = d/dt + u-V. In this form, one 
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recognizes the constitutive form of VCHE as a variant of the rate-dependent incompressible 
homogeneous fluid of second grade |15[ , [[Uj , whose viscous dissipation, however, is modified 
by the Helmholtz operator (1 — a 2 A). There is a tradition at least since Rivlin |I7j of mod- 



eling turbulence by using continuum mechanics principles such as objectivity and material 
frame indifference (see also [0). For example, this sort of approach is taken in deriving 



Reynolds stress algebraic equation models |19|| . Rate-dependent closure models of mean 



turbulence such as the VCHE have also been obtained by a two-scale DIA approach |2(| 
and by renormalization group methods |21] . 



I. Comparison of VCHE with LES and RANS models 

Reynolds-averaged Navier-Stokes (RANS) models of turbulence are part of the classic 
theoretical development of the subject p2|, pj, [Ell. The related Large Eddy Simulation 



(LES) turbulence modeling approach |[25|, ||26||, |27||, provides an operational definition of 



the intuitive idea of Eulerian resolved scales of motion in turbulent flow. In this approach 
a filtering function T(r) is introduced and the Eulerian velocity field U^ is filtered in an 
integral sense, as 

u(r)= / dV F{r - r') U B (r') . (33) 

JRZ 

This convolution of Ug with T defines the large scale, resolved, or filtered velocity, u. The 
corresponding small scale, or subgrid scale velocity, u', is then defined as the difference, 

u'(r) = U B (r)-u(r). (34) 

When this filtering operation is applied to the Navier-Stokes system, the following dynamical 
equation is obtained for the filtered velocity, u, cf. Eq. (|32|), 



d — 

— u + u-Vu = -divT- Vp+uAu, V-u = 0, (35) 

ot 
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in which p is the filtered pressure field (required to maintain V • u = 0) and the tensor 
difference 



T = (U s U £ )-uu, (36) 

represents the subgrid scale stress due to the turbulent eddies. This subgrid scale stress 
tensor appears in the same form as the Reynolds stress tensor obtained from Reynolds 
averaging the Navier-Stokes equation. 

The results of Chen et al. ||- |§, may be given either an LES, or RANS interpretation 



simply by comparing the constitutive form of the VCHE in (|32]) term by term with equation 
(J35|), provided one may ignore the difference between Eulerian mean, and Lagrangian mean 
velocities as being of higher order. Additional LES interpretations, discussions and numerical 
results for forced-turbulence simulations of the VCHE model, or the Navier-Stokes alpha 
model, will be presented below. 

III. DIRECT NUMERICAL SIMULATION, RELATED DEFINITIONS AND 

ENERGY SPECTRA 



In Fourier space, the viscous alpha model equation (|25f) with isotropic viscosity can be 
written for the Lagrangian mean velocity v as follows: 



5v k _ £,i,w„ „ ,.,J 



dt 



P(k)(uxq) k -z/£; 2 v k + f k , (37) 



k • v k = 0, (38) 

where P^ = Sij — kikj/k 2 is the incompressible projection operator. The Eulerian mean 
velocity satisfies u = (1 — a 2 A) -1 v (or u(k) = (1 + a 2 A; 2 ) -1 v(k)), q = V x v and f k is 
the forcing term for the k-th velocity component. We emphasize that 1/a acts as the cutoff 
wavenumber for the nonlinearity in alpha model. 
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A pseudo-spectral code has been developed for solving the above equations in a cubic 
box for periodic boundary conditions. The code is written in message-passing language for 
the SGI-ORIGIN-2000 machine in the Advanced Computing Laboratory at the Los Alamos 
National Laboratory The viscous term is integrated analytically in time. The other terms 
are discretized using a second-order Adams-Bashforth scheme. The nonlinear terms are 



calculated using pseudo-spectral methods j28fl . The time evolution of Eq. (j37|) can be 
written: 

-± ^P '- = P(k)[-(u x q)lexp(-vk 2 At) 

~ ^(v X q)r 1 exp(-2z/A; 2 At)] + / k . (39) 

In this paper, we adopt the following definitions: the mean velocity fluctuation, u', is 
defined as 



2 2 f°° 

u'=^J-E t = (-J E(k)dk) 



where E t is the total energy defined in Eq. (|27|) and E(k) is the energy spectrum: 

k+l/2 

E(k)= ^u(k')-v(k'). 

k-l/2 

The Taylor microscale and the mean dissipation rate are defined respectively as follows: 

A = (15z//e) 1/2 n', 

^■oo 

e = 2u k 2 E(k)dk. 

Jo 

The large eddy turn-over time is given by: 

v! 

where Lf is the integral length. The Kolmogorov dissipation scale r\ is defined as (u 3 /e)i, 

with corresponding wave number kj = 1/rj. The Taylor microscale Reynolds number is 

defined by 

u'X 

tt\ — . 

v 
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To maintain a statistical steady state, the forcing term fk was introduced in the first two 
shells of the Fourier modes (k < 2.5), in which the kinetic energy of each mode in the two 
shells was forced to be constant in time, while the energy ratio between shells consistent with 
£,-5/3 j n orc [ er t approximate a larger inertial sub-range ||28|| . Most simulations were carried 
out for about ten large eddy turnover times before recording any data. The initial velocity 
was a Gaussian field with a prescribed energy spectrum: ~ fc 4 exp[— (k/ko) 2 }, where ko is 
a constant whose value is approximately 5. Statistics were obtained by averaging physical 
quantities over several eddy turnover times. 

DNS has been carried out for three cases: with a = (the Navier-Stokes equations), 
1/32 and 1/8 for the same viscosity v = 0.001. The corresponding R\ are 147,182 and 
279, respectively. In Fig. 1, we show the energy spectra of the three simulations. We note 
that because each simulation used the same viscosity, the higher Reynolds number flows 
correspond to smaller Taylor microscales, leading to more compact energy spectra. This is 



similar to the results of subgrid simulations for higher Reynolds number flows ||29| . In Fig. 2, 
we compare the energy spectrum for a = 1/8 with mesh sizes of 256 3 and 64 3 . It is seen that 
the energy spectrum at the large scales (in the inertial range) are the same, indicating that 
the large scale flow properties for a = 1/8 can be preserved in a less-resolved simulation. 

To see how the energy cascades from the large scale modes to the small scale modes, in 
Fig. 3 we show the energy transfer spectrum, 

n(k)=v(k)-P(k)(uxq) k , 

as a function of the wave number, k. We note that the energy transfer spectrum for all 
three cases with different a's are constant and agree quite well in the inertial range for 
k < 20, except for the effect of noise. This result indicates that the alpha model preserves 



the fundamental properties of the Kolmogorov energy cascade in the inertial range |29 
For the dissipation range (k > 20), however, the change of the energy transfer spectrum is 
significant for a = 1/8. 
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IV. VORTICITY STRUCTURES AND ALIGNMENT 

As contrasted with the hyperviscosity approach 0,^1|] in which the normal dissipation 
operator is replaced by a higher order Laplace operator and therefore only small scale fluid 
motions are affected, in the alpha model the nonlinear vortex stretching dynamics is modified 
as shown in Eq. (p6|), where the vorticity q is defined by the curl of the velocity v while 
the advective velocity is u = (1 — a 2 A) _1 v. For a not zero, using a smoothed velocity u 
rather than the original velocity v suppresses the vortex stretching dynamics, especially for 
the small scale vortices whose size is less than a. 

To give a qualitative idea about how the vortex structures change with increasing a, in 
Fig. 4 (a-c) we present the iso-surfaces of vorticity for q/q' = 2 when a = 0, 1/32 and 1/8, 
where q' = (q 2 ) 1 ^ 2 is the root-mean- square value of q. It is seen that the tube-like vortex 
structures |32| persist in all three simulations, implying that the alpha model does not change 
the qualitative feature of stretching physics. On the other hand, it is evident that with 
increasing of a, the vortex aspect ratio (characteristic vortex radius/characteristic vortex 
length) decreases. A similar phenomenon has been noticed in hyperviscosity simulations 
also |30],[31j (where the vortex stretching dynamics is not suppressed) and in subgrid model 
simulations |2i| . However, the physical mechanisms for this phenomena are rather different 
in the alpha model from the actions of hyperviscosity, since the alpha term affects the 
nonlinearity (the cause of vortex stretching), while hyperviscosity does not. 

One of the most important properties of the vortex stretching dynamics in the Navier- 
Stokes turbulence is the so-called alignment phenomenon |53]], in which the three-dimensional 



vorticity is locally preferentially aligned with the direction of the second eigenvector of the 
symmetric strain-rate tensor, Sij = l/2(dui/dxj + duj/dxi). In three-dimensional space, Sij 
is a 3 x 3 matrix and has three eigenvalues, Ai, A 2 and A 3 . The incompressibility condition 
leads to: 

Ai + A 2 + A 3 = 0. 

Assume Ai > A 2 > A3, then Ai > and A3 < 0. For the Navier-Stokes turbulence, it has been 
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found |kj| that the spatially averaged A2 is greater than zero. Therefore in Navier-Stokes 
turbulence, the stretching dynamics of vorticity is dominated by two directional expansions 
with positive eigenvalues and one directional contraction with negative eigenvalue. Since 
the alpha model is equivalent to filtering the transport velocity in the stretching term for 
the vorticity equation Eq. Q26D, it is interesting to study how the alignment phenomenon 
changes in the alpha model when compared with the Navier-Stokes model. 

In Fig. 5 (a-c), we present the probability density functions (PDFs) of the cosine of 
the angles between the local vorticity and the eigenvectors of the local strain-rate tensor for 
a = 0, 1/32 and 1/8. The solid lines, the dotted lines and the dotted-dash lines correspond to 
the maximum eigenvalue, the middle eigenvalue and the minimum eigenvalue, respectively. 
As expected, the results in Fig. 5(a) (the Navier-Stokes case) are very similar to those results 
in [f33|1 . The major feature in this plot is that the PDF corresponding to the middle eigenvalue 
cos(#2) is peaked when cos (#2) = ±1? implying the alignment mentioned earlier. The PDF 
of cos(9 3 ) is peaked at the origin, indicating that the vorticity is locally perpendicular to 
the direction associated with the minimum eigenvalue. The PDF of cos{9\) is almost flat, 
implying that the direction of vorticity is essentially decorrelated from the direction of the 
maximum eigenvalue direction. With increasing a, the property of the minimum eigenvalue 
direction is qualitatively the same, but the PDF of cos{9\) starts forming a peak at ±1. 
For a = 1/8, the PDF values for cos(0i) at ±1 are even bigger than those of cos (#2). This 
result is quite different from the case of the Navier-Stokes equation. We suspect that this 
new alignment phenomenon (that the vorticity aligns with the direction of the maximum 
eigenvalue) is connected with the observation (as shown in Fig 4), that the high amplitude 
vorticity gets thicker as a increases, a result that has also been observed in other subgrid 
simulations [ f2"9f . 

To further quantify the change of eigenvalues as a function of a, in Fig. 6 (a-c), we show 
the PDFs of the eigenfunctions. There are two essential features in these plots: (i) as a 
increases, (A 2 ) keeps positivity, consistent with the result in the Navier-Stokes turbulence; 
(ii) all eigenvalues tend toward smaller values as a increases, meaning that the stretching is 
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being suppressed in the alpha model. To study the mean effect of the stretching term on the 
growth of the enstrophy, (q ■ q) , in Fig. 7 we present the PDF of Z = q ■ S ■ q as a function 
of a. The PDFs of these three quantities show very strong intermittency. The flatnesses are 
F qiSi . qj = 1965 for a = 0, 3566 for a = 1/32 and 5319 for a = 1/8. The PDFs are strongly 
positively skewed (24.28,35.71 and 47.02) which is consistent with the expectation that the 
stretching term in the vorticity equation gives a net positive contribution to the dynamics 
of enstrophy. 

V. TWO-POINT STATISTICS 

Two-point statistics, primarily the scaling relations of the velocity structure functions 
in the inertial range, were proposed by Kolmogorov in 1941 |34| based on the self-similarity 



hypothesis and in 1962 |35| using the idea of the refined similarity hypothesis. Understanding 
the inertial range intermittent dynamics of fluid turbulence has been a very important 
focal point for the last few decades, and the two-point inertial range statistics have been 
extensively examined in experiments at high Reynolds numbers and by numerical simulations 
at low to moderate Reynolds numbers [ 31] . In this section, we study the probability density 



function of the velocity increment and the scalings of the velocity structure functions for 
systems with various values of a using DNS data, and then compare our measurements with 
the Kolmogorov theories. 

In Fig. 8, we compare the PDFs of the longitudinal velocity gradient, du/dx, for various 
alpha values. The derivative here is calculated by a central difference in physical space and a 
spatial averaging is used for the ensemble averaging in the normalization factor calculation. 
The main information in this plot is that the PDFs for all three alpha values have close to 
an exponential tail for large amplitude events. However, the flatnesses of the PDFs (5.0, 
4.1 and 3.8 for a = 0, 1/32 and 1/8, respectively) decrease with increasing a. The observed 
reduction of small scale intermittency is consistent with the energy spectra in Fig. 1 where 
the high k energy spectra are truncated for large a values. The skewnesses of the PDFs 



are —0.48, —0.51 and —0.54. In Fig. 9, we present the normalized PDFs of the velocity 
increment , 

A r u = u(x + r) — u(x) 

in the inertial range when the separation r = 20 (in mesh units). It is seen that all three 
PDFs collapse quite well for most values of A r u, implying that the fundamental features of 
the two-point statistics in the inertial range do not change much as a varies. This is a very 
desirable feature in using the alpha model as a subgrid model. 

In Fig. 10, we show the second-order structure function, (A r u), versus r at different 
alpha values. It is seen that when r less than 5 mesh points, the Taylor expansion gives 
a simple scaling: (A r u) ~ r 2 . In the narrow inertial range (20 < r < 60), the simulation 



results with a = and 1/32 agree qualitatively with the Kolmogorov 2/3 scaling f3~4" |. The 
slightly slower than 2/3 growth may indicate the intermittent correction. For the case of 
a = 1/8, the filtering of small scale motions significantly decreases the structure function 
in the inertial range. For this case, the 2/3 scaling can only been seen for the very narrow 
region 40 < r < 60. 

In Fig. 11, we present the flatness of the velocity increment, ((A r u) 4 ) / langle(A r u) 2 ) 2 , 
versus the normalized separation, t/tq. The normalization length is taken to be half the 
box length tt. As expected, in the small scale region with r/r Q < 0.005, the alpha model 
is quite different from the Navier Stokes dynamics and the Navier-Stokes equation has the 
largest flatness value, implying that the Navier-Stokes turbulence is more intermittent than 
the alpha model dynamics, consistent with the vortex visualization in Fig. 4. It is very 
interesting to note, however, that the flatnesses versus r /r$ are almost identical for all three 
alpha values studied in the narrow inertial range. The decaying of the flatness as a function 
of r/r implies the existence of an intermittency correction for the scaling exponents. The 
numerical measurement gives ((A r u) 4 ) / ((A r u) 2 ) 2 ~ r -° 12 j n the inertial range. This scaling 
exponent agrees qualitatively with the value of —0.11 from other direct numerical simulation 
results |3l| . 
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VI. CONCLUDING REMARKS 

Our main objective in this paper has been to investigate the utility of the recently 
proposed alpha equations as a subgrid scale model for three- dimensional turbulence. Our 
main tool has been direct numerical simulation (DNS) of turbulence in a periodic box. Our 
conclusions are based on comparisons of DNS using the Navier-Stokes alpha equations for 
several values of the alpha parameter, including the limiting case a = in which the Navier- 
Stokes equations are recovered. Our principal conclusion is that the alpha model simulations 
can reproduce most of the large scale features of Navier-Stokes turbulence even when these 
simulations do not resolve the fine scale dynamics, at least in the case of turbulence in a 
periodic box. 

We summarized known analytic properties of the alpha models, including an outline of 
their derivation and the associated assumptions, their simplification for the case of constant 
alpha, and their conservation properties. We also offered interpretations of nonlinear dy- 
namics of the alpha models and indicated the changes one might expect from the dynamics 
of the Navier-Stokes equations. 

One of our principal computational results is shown in Fig. 3, where two DNS simulations 
at a = 1/8, but using two different resolutions, are compared with each other and with 
a Navier-Stokes simulation (representing truth). The spectra are identical for the larger 
wavelengths, demonstrating (in this case) that one does not need to resolve the small scale 
dynamics to reproduce the large scale features of the turbulence. 

We further compared vorticity structures and alignment, and also two point statistics 
to illustrate the altered dynamics of the alpha models. In general, we found consistency of 
these DNS results with our expectations based on analysis. 

The use of the alpha equations as a subgrid model of turbulent flows in more complicated 
geometries and forcings remains to be studied. However we believe these first results are 
very promising. 
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VIII. FIGURE CAPTIONS 

Fig. 1. The energy spectrum, E(k), versus the wave number k for a = (solid line), 
1/32 (dotted line) and 1/8 (dotted-dash line). In the inertial range (k < 10), a power 
spectrum with A;" 5 / 3 can be identified. 

Fig. 2. Energy spectra for a — 1/8 with 256 3 (dotted line) and 64 3 (dashed-dot line), and 
compared with the Navier-Stokes simulation (the solid line). 

Fig. 3. The energy transfer spectrum, Tl(k) versus k for a = (solid line), 1/32 (dotted 
line) and 1/8 (dotted-dash line). 

Fig. 4. Three-dimensional view of vorticity iso-surfaces for a = (a), 1/32 (b) and 1/8 (c) 
when q/q' = 2, where q' is the rms of the vorticity. The data are from 256 3 simulations 
and only one eighth of the simulation domain is shown. 

Fig. 5. The probability density function of cos(8) for simulations with a = (a), 1/32 (b) 
and 1/8 (c). In each plot, the solid line, the dotted line and the dotted-dash line are 
for the maximum, middle and minimum eigenvalue, respectively. 

Fig. 6. The probability density functions for the maximum eigenvalue Ai (a), the middle 
eigenvalue A2 (b) and the minimum eigenvalue A3 (c). The solid line, the dotted line 
and the dotted-dash line are for a — 0, 1/32 and 1/8, respectively. 

Fig. 7. Normalized PDFs of z — qiSijqj for a = (solid line), 1/32 (dotted line) and 1/8 
(dotted-dash line). 
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Fig. 8. Normalized PDFs of the longitudinal velocity gradient, du/dx for a = (solid line), 
1/32 (dotted line) and 1/8 (dotted-dash line). 

Fig. 9. Normalized PDFs of the velocity increment in the inertial range for a = (solid 
line), 1/32 (dotted line) and 1/8 (dotted-dash line). 

Fig. 10. The second order structure function as a function of r for a = (solid line), 1/32 
(dotted line) and 1/8 (dotted-dash line). The dashed line is the scaling prediction by 
Kolmogorov |34| . 



Fig. 11. The flatness, ((A r M) 4 )/((A r M) 2 ) 2 , as a function of separation r/r for a = (solid 
line), 1/32 (dotted line) and 1/8 (dotted-dash line). Here r is taken as half the box 
size, 7r. 
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